Chapter 17 - Exercises¶

In [1]:
library(tidyverse)
library(bayesrules)
library(bayesplot)
library(rstan)
library(rstanarm)
library(broom.mixed)
library(tidybayes)
library(forcats)
library(lme4)
# use multiple cores, otherwise hierarchical models with random slopes and intercepts are terribly slow
options(mc.cores=parallel::detectCores())

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.2     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.2     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
This is bayesplot version 1.10.0

- Online documentation and vignettes at mc-stan.org/bayesplot

- bayesplot theme set to bayesplot::theme_default()

   * Does _not_ affect other ggplot2 plots

   * See ?bayesplot_theme_set for details on theme setting

Loading required package: StanHeaders

rstan (Version 2.21.8, GitRev: 2e1f913d3ca3)

For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores()).
To avoid recompilation of unchanged Stan programs, we recommend calling
rstan_options(auto_write = TRUE)


Attaching package: ‘rstan’


The following object is masked from ‘package:tidyr’:

    extract


Loading required package: Rcpp

This is rstanarm version 2.21.4

- See https://mc-stan.org/rstanarm/articles/priors for changes to default priors!

- Default priors may change, so it's safest to specify priors, even if equivalent to the defaults.

- For execution on a local, multicore CPU with excess RAM we recommend calling

  options(mc.cores = parallel::detectCores())


Attaching package: ‘rstanarm’


The following object is masked from ‘package:rstan’:

    loo


Loading required package: Matrix


Attaching package: ‘Matrix’


The following objects are masked from ‘package:tidyr’:

    expand, pack, unpack

Exercise 17.6¶

In [2]:
head( sleepstudy )
A data.frame: 6 × 3
ReactionDaysSubject
<dbl><dbl><fct>
1249.56000308
2258.70471308
3250.80062308
4321.43983308
5356.85194308
6414.69015308
In [3]:
dim( sleepstudy )
  1. 180
  2. 3
In [4]:
sleepstudy  %>% summarize( subjects=n_distinct(Subject) )
A data.frame: 1 × 1
subjects
<int>
18

a)¶

We should group by the subject of the experiment. Since we expect that average reaction time differs from subject to subject, we should incorporate this information.

b)¶

See exercises 17.1 (b) and (a)

c)¶

Model 1 assumes that the slope of the decline (increase) of reaction time with the number of sleep-deprived days is the same for all subjects (by fixing a global slope). Model 2 does not make this assumption - the slope of the decline of reaction time is individual to the subject.

d)¶

In [5]:
options(repr.plot.width=15, repr.plot.height=5)
ggplot(sleepstudy, aes(x=Days, y=Reaction, group=Subject)) + 
    geom_smooth(method="lm", se=FALSE, linewidth=0.5)

`geom_smooth()` using formula = 'y ~ x'

The variance within the slopes for the individual slopes is significant, but not very large. In view of the very small amount of different subjects (18), model 1 might perform better.

Exercise 17.7¶

a)¶

Model 1:

In [6]:
sleep_model_1 <- stan_glmer(
  Reaction ~ Days + (1 | Subject), 
  data = sleepstudy, family = gaussian,
  prior_intercept = normal(250, 10, autoscale = TRUE),
  prior = normal(0, 1, autoscale = TRUE), 
  prior_aux = exponential(1, autoscale = TRUE),
  prior_covariance = decov(reg = 1, conc = 1, shape = 1, scale = 1),
  chains = 4, iter = 5000*2, seed = 84735
)

MCMC diagnostics:

In [7]:
options(repr.plot.width=15, repr.plot.height=10)
mcmc_trace(sleep_model_1)
mcmc_dens_overlay(sleep_model_1)
neff_ratio(sleep_model_1)
rhat(sleep_model_1)
(Intercept)
0.14215
Days
0.84515
b[(Intercept) Subject:308]
0.2416
b[(Intercept) Subject:309]
0.23105
b[(Intercept) Subject:310]
0.2283
b[(Intercept) Subject:330]
0.2191
b[(Intercept) Subject:331]
0.22415
b[(Intercept) Subject:332]
0.2195
b[(Intercept) Subject:333]
0.22085
b[(Intercept) Subject:334]
0.2156
b[(Intercept) Subject:335]
0.2367
b[(Intercept) Subject:337]
0.2293
b[(Intercept) Subject:349]
0.2203
b[(Intercept) Subject:350]
0.2301
b[(Intercept) Subject:351]
0.2225
b[(Intercept) Subject:352]
0.2282
b[(Intercept) Subject:369]
0.21445
b[(Intercept) Subject:370]
0.21545
b[(Intercept) Subject:371]
0.21495
b[(Intercept) Subject:372]
0.2184
sigma
0.6223
Sigma[Subject:(Intercept),(Intercept)]
0.20645
(Intercept)
1.00082664767411
Days
1.0000622494225
b[(Intercept) Subject:308]
1.00052999092531
b[(Intercept) Subject:309]
1.00083026408339
b[(Intercept) Subject:310]
1.00129926473946
b[(Intercept) Subject:330]
1.00056697528007
b[(Intercept) Subject:331]
1.00053645231849
b[(Intercept) Subject:332]
1.00095837241386
b[(Intercept) Subject:333]
1.00039350335963
b[(Intercept) Subject:334]
1.00032112854831
b[(Intercept) Subject:335]
1.00070866485161
b[(Intercept) Subject:337]
1.00038691166775
b[(Intercept) Subject:349]
1.00088752697069
b[(Intercept) Subject:350]
1.00053319383892
b[(Intercept) Subject:351]
1.0005064631744
b[(Intercept) Subject:352]
1.00072584724475
b[(Intercept) Subject:369]
1.00097258614422
b[(Intercept) Subject:370]
1.00064792950943
b[(Intercept) Subject:371]
1.0006000997014
b[(Intercept) Subject:372]
1.00071089475644
sigma
0.999899120713781
Sigma[Subject:(Intercept),(Intercept)]
1.00160149418865

Looks good!

Model 2:

In [8]:
sleep_model_2 <- stan_glmer(
  Reaction ~ Days + (Days | Subject), 
  data = sleepstudy, family = gaussian,
  prior_intercept = normal(250, 10, autoscale = TRUE),
  prior = normal(0, 1, autoscale = TRUE), 
  prior_aux = exponential(1, autoscale = TRUE),
  prior_covariance = decov(reg = 1, conc = 1, shape = 1, scale = 1),
  chains = 4, iter = 5000*2, seed = 84735
)

MCMC diagnostics:

In [9]:
options(repr.plot.width=15, repr.plot.height=10)
mcmc_trace(sleep_model_2)
mcmc_dens_overlay(sleep_model_2)
neff_ratio(sleep_model_2)
rhat(sleep_model_2)
(Intercept)
0.55455
Days
0.36975
b[(Intercept) Subject:308]
0.7321
b[Days Subject:308]
0.55325
b[(Intercept) Subject:309]
0.98435
b[Days Subject:309]
0.7357
b[(Intercept) Subject:310]
0.92715
b[Days Subject:310]
0.6761
b[(Intercept) Subject:330]
0.698
b[Days Subject:330]
0.56885
b[(Intercept) Subject:331]
0.7244
b[Days Subject:331]
0.5774
b[(Intercept) Subject:332]
0.91605
b[Days Subject:332]
0.6542
b[(Intercept) Subject:333]
0.90815
b[Days Subject:333]
0.65305
b[(Intercept) Subject:334]
0.87285
b[Days Subject:334]
0.6307
b[(Intercept) Subject:335]
0.6265
b[Days Subject:335]
0.53725
b[(Intercept) Subject:337]
0.91165
b[Days Subject:337]
0.7212
b[(Intercept) Subject:349]
0.7835
b[Days Subject:349]
0.5859
b[(Intercept) Subject:350]
0.68095
b[Days Subject:350]
0.523
b[(Intercept) Subject:351]
0.8475
b[Days Subject:351]
0.6517
b[(Intercept) Subject:352]
0.9854
b[Days Subject:352]
0.6638
b[(Intercept) Subject:369]
0.8849
b[Days Subject:369]
0.6092
b[(Intercept) Subject:370]
0.6637
b[Days Subject:370]
0.54355
b[(Intercept) Subject:371]
0.94915
b[Days Subject:371]
0.6303
b[(Intercept) Subject:372]
1.00145
b[Days Subject:372]
0.62705
sigma
0.77415
Sigma[Subject:(Intercept),(Intercept)]
0.4602
Sigma[Subject:Days,(Intercept)]
0.26435
Sigma[Subject:Days,Days]
0.4476
(Intercept)
1.00021233355127
Days
1.00012363141835
b[(Intercept) Subject:308]
0.999964511619639
b[Days Subject:308]
1.0000334938779
b[(Intercept) Subject:309]
0.999960299235084
b[Days Subject:309]
0.999905048719386
b[(Intercept) Subject:310]
1.00004840751115
b[Days Subject:310]
0.999933151636817
b[(Intercept) Subject:330]
1.00022082980943
b[Days Subject:330]
1.00020220817309
b[(Intercept) Subject:331]
1.00002093496713
b[Days Subject:331]
1.00009787002183
b[(Intercept) Subject:332]
0.999969137141594
b[Days Subject:332]
1.00015545536961
b[(Intercept) Subject:333]
1.00003011196354
b[Days Subject:333]
0.999984943519625
b[(Intercept) Subject:334]
1.00013628184625
b[Days Subject:334]
1.00006951334323
b[(Intercept) Subject:335]
0.999943627457239
b[Days Subject:335]
1.0000816868977
b[(Intercept) Subject:337]
1.00011821131814
b[Days Subject:337]
0.999890193456099
b[(Intercept) Subject:349]
1.0001444437165
b[Days Subject:349]
1.00010805015095
b[(Intercept) Subject:350]
1.00054919237587
b[Days Subject:350]
1.00018046907854
b[(Intercept) Subject:351]
1.00004289114336
b[Days Subject:351]
1.00014495586934
b[(Intercept) Subject:352]
0.999865177134585
b[Days Subject:352]
0.999975394078911
b[(Intercept) Subject:369]
1.0002859705965
b[Days Subject:369]
0.999958326829176
b[(Intercept) Subject:370]
1.00008253180516
b[Days Subject:370]
1.00001082175391
b[(Intercept) Subject:371]
0.999979469653042
b[Days Subject:371]
1.00006240655609
b[(Intercept) Subject:372]
1.00003800077104
b[Days Subject:372]
1.00003196589408
sigma
1.00029493215927
Sigma[Subject:(Intercept),(Intercept)]
1.00034192996102
Sigma[Subject:Days,(Intercept)]
1.00018938007471
Sigma[Subject:Days,Days]
1.00027339003668

Looks good!

b)¶

In [10]:
tidy(sleep_model_2, effects = "fixed", conf.int = TRUE, conf.level = 0.80)
A tibble: 2 × 5
termestimatestd.errorconf.lowconf.high
<chr><dbl><dbl><dbl><dbl>
(Intercept)251.437526.560464242.769578259.98264
Days 10.405321.687297 8.167352 12.55381

Posterior median regression model:

Reaction in $\mu$s = 251.4 + 10.4 * Days

c)¶

Get posterior samples:

In [11]:
set.seed(84735)
sleep_chains_2 <- sleep_model_2 %>%
    spread_draws(`(Intercept)`, b[term, Subject], `Days`) %>% 
    pivot_wider(names_from = term, names_glue = "b_{term}", values_from = b) %>% 
    mutate(subject_intercept = `(Intercept)` + `b_(Intercept)`, subject_Days = Days+b_Days)

head( sleep_chains_2 )
A grouped_df: 6 × 10
.chain.iteration.draw(Intercept)SubjectDaysb_(Intercept)b_Dayssubject_interceptsubject_Days
<int><int><int><dbl><chr><dbl><dbl><dbl><dbl><dbl>
111251.3246Subject:3089.111127 10.30740 9.7072159261.632018.818343
111251.3246Subject:3099.111127-32.91128-7.4075791218.4133 1.703548
111251.3246Subject:3109.111127-50.59694-1.1558622200.7276 7.955265
111251.3246Subject:3309.111127 30.11841-2.6720126281.4430 6.439114
111251.3246Subject:3319.111127 15.74251 0.5656770267.0671 9.676804
111251.3246Subject:3329.111127 19.96727-0.7068745271.2919 8.404252

80% credible intervals for Days regression coefficient for each subject:

In [12]:
sleep_chains_2  %>% 
    group_by( Subject ) %>% 
    summarize( Days_lower=quantile(subject_Days, 0.1), Days_upper=quantile(subject_Days, 0.9) ) %>% 
    arrange( Subject )
A tibble: 18 × 3
SubjectDays_lowerDays_upper
<chr><dbl><dbl>
Subject:30816.52413522.819969
Subject:309-1.911605 4.418053
Subject:310 1.242503 7.654293
Subject:330 2.750641 9.210461
Subject:331 4.52310610.821231
Subject:332 7.36001313.257076
Subject:333 7.42404913.463109
Subject:334 8.49560914.362475
Subject:335-3.498164 2.962450
Subject:33716.40286122.687440
Subject:349 8.11300314.405007
Subject:35013.64967720.083856
Subject:351 4.56402910.438934
Subject:35211.31243117.221914
Subject:369 8.40093614.273878
Subject:37011.65600818.218000
Subject:371 6.57480612.373092
Subject:372 9.00654514.792631

d)¶

In [13]:
tidy(sleep_model_2, effects = "ran_pars")
A tibble: 4 × 3
termgroupestimate
<chr><chr><dbl>
sd_(Intercept).Subject Subject 24.20655929
sd_Days.Subject Subject 6.91446511
cor_(Intercept).Days.SubjectSubject 0.08132467
sd_Observation.Residual Residual25.97541320

$\sigma_y \approx 26.0$, $\sigma_0 \approx 24.2$, $\sigma_1 \approx 6.9$, $\rho \approx 0.08$

Per measurement, the reaction times vary with a standard deviation of 26 $\mu$s. The model intercept ($\sigma_0$) explains much more of the variance (93.4%) than the slope coefficient:

In [14]:
pi0 = 26^2/(26^2+6.9^2)
pi0
0.93420488937411
In [15]:
pi1 = 6.9^2/(26^2+6.9^2)
pi1
0.0657951106258896

The correlation between intercept and slope is moderatly positive, indicating that people with higher baseline reaction times respond more strongly to sleep deprivation.

Exercise 17.8¶

a)¶

Posterior median regression models:

In [16]:
sleep_chains_2  %>% 
    group_by( Subject ) %>% 
    summarize( beta0=median(subject_intercept), beta1=median(subject_Days) ) %>% 
    arrange( beta1 )
A tibble: 18 × 3
Subjectbeta0beta1
<chr><dbl><dbl>
Subject:335250.3482-0.2116396
Subject:309214.6996 1.2428889
Subject:310216.1619 4.4109275
Subject:330272.3174 6.0482060
Subject:351255.3197 7.5636197
Subject:331271.1705 7.7641621
Subject:371252.0949 9.4645675
Subject:332259.303010.3602429
Subject:333266.458510.5063303
Subject:349229.053311.1931523
Subject:369254.458911.3539589
Subject:334245.280011.3734917
Subject:372262.596311.9242890
Subject:352270.489014.2736235
Subject:370228.869014.8111600
Subject:350240.255916.7603131
Subject:337283.290819.5673004
Subject:308254.305919.5936416

It appears that Subject 335 even improved slightly with sleep deprivation:

Reaction = 250.3 - 0.21 * Days

In [17]:
options(repr.plot.width=15, repr.plot.height=5)
ggplot( sleepstudy %>% filter(Subject==335), aes(x=Days, y=Reaction) ) + 
    geom_line() + 
    geom_point() +
    geom_smooth( method="lm" )

`geom_smooth()` using formula = 'y ~ x'

(notice the shrinkage)

b)¶

Subject 308: Reaction =254.3 + 19.6*Days

In [18]:
options(repr.plot.width=15, repr.plot.height=5)
ggplot( sleepstudy %>% filter(Subject==308), aes(x=Days, y=Reaction) ) + 
    geom_line() + 
    geom_point() +
    geom_smooth( method="lm" )

`geom_smooth()` using formula = 'y ~ x'

c)¶

In [19]:
sleep_chains_2  %>% 
    group_by( Subject ) %>% 
    summarize( beta0=median(subject_intercept), beta1=median(subject_Days) ) %>% 
    arrange( beta0 )
A tibble: 18 × 3
Subjectbeta0beta1
<chr><dbl><dbl>
Subject:309214.6996 1.2428889
Subject:310216.1619 4.4109275
Subject:370228.869014.8111600
Subject:349229.053311.1931523
Subject:350240.255916.7603131
Subject:334245.280011.3734917
Subject:335250.3482-0.2116396
Subject:371252.0949 9.4645675
Subject:308254.305919.5936416
Subject:369254.458911.3539589
Subject:351255.3197 7.5636197
Subject:332259.303010.3602429
Subject:372262.596311.9242890
Subject:333266.458510.5063303
Subject:352270.489014.2736235
Subject:331271.1705 7.7641621
Subject:330272.3174 6.0482060
Subject:337283.290819.5673004

Subject 337: Reaction = 283.3 + 19.6 * Days

In [20]:
options(repr.plot.width=15, repr.plot.height=5)
ggplot( sleepstudy %>% filter(Subject==337), aes(x=Days, y=Reaction) ) + 
    geom_line() + 
    geom_point() +
    geom_smooth( method="lm" )

`geom_smooth()` using formula = 'y ~ x'

d)¶

Subject 309: Reaction = 214.7 + 1.24 * Days

In [21]:
options(repr.plot.width=15, repr.plot.height=5)
ggplot( sleepstudy %>% filter(Subject==309), aes(x=Days, y=Reaction) ) + 
    geom_line() + 
    geom_point() +
    geom_smooth( method="lm" )

`geom_smooth()` using formula = 'y ~ x'

e)¶

Without posterior_predict() for subject 308:

In [22]:
set.seed(84735)
sleep_model_2_df <- data.frame( sleep_model_2 ) 
sleep_model_2_df %>% 
    transmute( 
        subject_intercept = X.Intercept.+b..Intercept..Subject.308.,
        subject_slope = Days+b.Days.Subject.308.,
        sigma_y = sigma
    ) %>% 
    mutate( mu = subject_intercept+5*subject_slope ) %>% 
    mutate( ynew = rnorm(nrow(sleep_model_2_df), mean=mu, sd=sigma_y) ) %>% 
    summarize( ynew_median=median(ynew), ynew_lower=quantile(ynew, 0.1), ynew_upper=quantile(ynew, 0.9))
A data.frame: 1 × 3
ynew_medianynew_lowerynew_upper
<dbl><dbl><dbl>
351.7238317.564387.0963

Doing this for a new subject is a bit more involved (extracting covariance matrix coefficients and sampling from multivariate normal distribution), I leave this up to you.

With posterior_predict():

In [23]:
predict_reaction_5days <- posterior_predict(
  sleep_model_2, 
  newdata = data.frame(Subject = c("308", "myself"), Days = c(5,5))
)

Median values:

In [24]:
data.frame( predict_reaction_5days ) %>% 
    rename( subject_308=X1, myself=X2 ) %>% 
    summarize_all( median )
A data.frame: 1 × 2
subject_308myself
<dbl><dbl>
352.1544303.5474

Distributions:

In [25]:
mcmc_areas(predict_reaction_5days, prob = 0.8) +
 ggplot2::scale_y_discrete(labels = c("Subject 308", "myself"))

Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

Exercise 17.9¶

a)¶

Are they wrong?

In [26]:
pp_check(sleep_model_1) +
  xlab("Days")
pp_check(sleep_model_2) +
  xlab("Days")

Model 2 looks slightly better, but not much! Both are not perfect, not really capturing the skewness in the data.

Are they fair? Hopefully the subjects were not harmed in long-term by the study - they probably signed some sort of waiver anyway. I can't imagine how the results of the study could be used to discriminate somebody, the results are very general. So: The model appears to be fair.

How accurate are the posterior predictions?

In [27]:
prediction_summary(model = sleep_model_1, data = sleepstudy)
A data.frame: 1 × 4
maemae_scaledwithin_50within_95
<dbl><dbl><dbl><dbl>
16.727880.50883830.62222220.9611111
In [28]:
prediction_summary(model = sleep_model_2, data = sleepstudy)
A data.frame: 1 × 4
maemae_scaledwithin_50within_95
<dbl><dbl><dbl><dbl>
12.417650.44415260.70.9777778

Model 2 appears to be a bit better in median absolute error.

b)¶

In [29]:
elpd1 <- loo(sleep_model_1, k_threshold=0.7)
elpd2 <- loo(sleep_model_2, k_threshold=0.7)
loo_compare(elpd1, elpd2)

All pareto_k estimates below user-specified threshold of 0.7. 
Returning loo object.

3 problematic observation(s) found.
Model will be refit 3 times.


Fitting model 1 out of 3 (leaving out observation 8)


Fitting model 2 out of 3 (leaving out observation 57)


Fitting model 3 out of 3 (leaving out observation 60)
A compare.loo: 2 × 8 of type dbl
elpd_diffse_diffelpd_loose_elpd_loop_loose_p_loolooicse_looic
sleep_model_2 0.00000 0.00000-862.774422.9506834.853789.0869651725.54945.90137
sleep_model_1-22.3128512.18954-885.087314.4016819.547773.3862031770.17528.80337

-22.3 $\pm$ 2 * 12.2 = (-46.7, 2.1). An ELPD difference of zero is within two standard errors of the estimated difference. There is no strong evidence that model 2 is better than model 2, thus we should favour the less flexible model 1 (also in light of the limited amount of data).

Thus I prefer a hierarchical model with random intercepts, but slopes, i.e. I assume that all subjects react in the same way to sleep deprivation.

Exercise 17.10¶

a)¶

\begin{eqnarray*} Y_{ij}|\beta_{0j},\beta_1,\sigma_y &\sim& N(\mu_{ij}, \sigma_y^2) \quad \text{with} \quad \mu_{ij} = \beta_{0j} + \beta_{1} X _{ij}\\ \beta_{0j} |\beta_0, \sigma_0 &\sim & N\left( \beta_0, \sigma_0^2 \right)\\ \beta_{0c} &=& N(m_0, s_0^2) \\ \beta_1 &=& N(0, s_1^2)\\ \sigma_y &=& \text{Exp}(l_y) \\ \sigma_0 &=& \text{Exp}(l_0) \end{eqnarray*}

b)¶

To interpret the coefficients, we have to make an assumption about the encoding of the two levels of $X_{ij}$. Let me assume that polite is the baseline. Then $\beta_{0j}$ is the mean polite pitch of subject $j$ and $\beta_0$ is the mean polite pitch among all the subjects.

c)¶

$\sigma_y$ is the within-group variance, i.e. the unexplained spread of polite pitch for one subject that is still remaining after fitting the linear model. $\sigma_0$ is the between-group variance, i.e. how much polite pitch varies between different subjects.

Exercise 17.11¶

a)¶

In [30]:
head( voices )
A tibble: 6 × 4
subjectscenarioattitudepitch
<fct><fct><fct><dbl>
AApolite 229.7
AAinformal237.3
ABpolite 236.8
ABinformal251.0
ACpolite 267.0
ACinformal266.0
In [31]:
dim( voices )
  1. 84
  2. 4

Subjects:

In [32]:
voices %>% summarize( subjects=n_distinct(subject) )
A tibble: 1 × 1
subjects
<int>
6

This is a rather small number! It definitely makes sense to choose the simpler model.

Dialogs per subject:

In [33]:
voices %>% na.omit() %>% count( subject )
A tibble: 6 × 2
subjectn
<fct><int>
A14
B13
C14
D14
E14
F14

14 dialogs per subject (with one exception).

b)¶

In [34]:
ggplot( voices %>% na.omit() ) + geom_boxplot( aes(x=subject, y=pitch, color=attitude) )

The between-group variance appears to be larger than the within-group variance (what makes sense, my pitch will not change so much that I sound like somebody else).

Exercise 17.12¶

a)¶

Simulate the model:

In [35]:
voice_model <- stan_glmer(
  pitch ~ attitude + (1 | subject), 
  data = voices, family = gaussian,
  prior_intercept = normal(250, 10, autoscale = TRUE),
  prior = normal(0, 1, autoscale = TRUE), 
  prior_aux = exponential(1, autoscale = TRUE),
  prior_covariance = decov(reg = 1, conc = 1, shape = 1, scale = 1),
  chains = 4, iter = 5000*2, seed = 84735
)

Trace plot:

In [36]:
options(repr.plot.width=15, repr.plot.height=10)
mcmc_trace(voice_model)
mcmc_dens_overlay(voice_model)
mcmc_acf(voice_model)
neff_ratio(voice_model)
rhat(voice_model)

Warning message:
“The `facets` argument of `facet_grid()` is deprecated as of ggplot2 2.2.0.
ℹ Please use the `rows` argument instead.
ℹ The deprecated feature was likely used in the bayesplot package.
  Please report the issue at <https://github.com/stan-dev/bayesplot/issues/>.”
(Intercept)
0.2628
attitudepolite
0.8317
b[(Intercept) subject:A]
0.2719
b[(Intercept) subject:B]
0.2759
b[(Intercept) subject:C]
0.27375
b[(Intercept) subject:D]
0.27585
b[(Intercept) subject:E]
0.27395
b[(Intercept) subject:F]
0.27165
sigma
0.6727
Sigma[subject:(Intercept),(Intercept)]
0.3645
(Intercept)
1.00056961772388
attitudepolite
0.999903917022183
b[(Intercept) subject:A]
1.0005005751095
b[(Intercept) subject:B]
1.00048034228008
b[(Intercept) subject:C]
1.00055133086772
b[(Intercept) subject:D]
1.000517849664
b[(Intercept) subject:E]
1.00047154616439
b[(Intercept) subject:F]
1.00056747020481
sigma
1.00002219680062
Sigma[subject:(Intercept),(Intercept)]
1.00121714449572

Looks good!

Posterior predictive check:

In [37]:
options(repr.plot.width=15, repr.plot.height=5)
pp_check(voice_model) + xlab("Attitude")

A lot of variance, but looks ok, no fundamental systematic error.

b)¶

In [38]:
tidy(voice_model, effects = "fixed", conf.int = TRUE, conf.level = 0.95)
A tibble: 2 × 5
termestimatestd.errorconf.lowconf.high
<chr><dbl><dbl><dbl><dbl>
(Intercept) 203.1464924.600019148.61465256.131177
attitudepolite-19.31093 6.427638-32.14464 -6.691067

(notice that the baseline was chosen different than in the theoretical explanations above)

$\beta_0 \in [148.6,256.1]$

c)¶

$\beta_1 \in [-32.1,-6.7]$

d)¶

The 95% credible interval for $\beta_1$ is quite large, but does not contain 0, pitch appears therefore to be significantly correlated with attitude.

Exercise 17.13¶

a)¶

pitch = 203.1 - 19.31 * $X$

where $X = \left\{\begin{array}{cc} 0, & \text{informal}\\ 1, & \text{polite} \end{array} \right.$

b)¶

In [39]:
subject_summary <- voice_model %>%
  spread_draws(`(Intercept)`, b[,subject]) %>% # extract dataframe
  mutate(subject_intercept = `(Intercept)` + b) %>% # compute individual runner intercepts
  select(-`(Intercept)`, -b) %>% # drop global intercept and b
  median_qi(.width = 0.80) %>%  # compute credible intervals
  select(subject, subject_intercept, .lower, .upper)

subject_summary
A tibble: 6 × 4
subjectsubject_intercept.lower.upper
<chr><dbl><dbl><dbl>
subject:A259.3522248.3940270.2947
subject:B155.9381144.8589167.0996
subject:C240.9517230.1386252.1316
subject:D113.6184102.7137124.8820
subject:E266.5357255.6265277.4952
subject:F179.0661168.1452189.9411

subject A: pitch = 259.4 - 19.31 * $X$

subject B: pitch = 155.9 - 19.31 * $X$

c)¶

In [40]:
set.seed(84735)
predict_new_scenario <- posterior_predict(
  voice_model, 
  newdata = data.frame(subject = c("A", "F", "myself"), attitude = c("polite","polite", "polite"))
)

Median values:

In [41]:
data.frame( predict_new_scenario ) %>% 
    rename( A=X1, B=X2, myself=X3 ) %>% 
    summarize_all( median )
A data.frame: 1 × 3
ABmyself
<dbl><dbl><dbl>
240.0855159.5334183.9078

Distributions:

In [42]:
options(repr.plot.width=15, repr.plot.height=5)
mcmc_areas(predict_new_scenario, prob = 0.8) +
 ggplot2::scale_y_discrete(labels = c("A", "F", "myself"))

Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

Since the between-group variability is larger than the within-group variability, as soon as the baseline pitch is known, the distributions are quite narrow (A and F). If the pitch is not known, the distributions are much wider (myself).

Exercise 17.14¶

Exploratory analysis¶

In [43]:
sample_n( coffee_ratings_small, 10 )
A tibble: 10 × 11
farm_nametotal_cup_pointsaromaflavoraftertasteaciditybodybalanceuniformitysweetnessmoisture
<fct><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl>
NA 82.587.587.507.427.587.677.4210100.11
NA 82.337.337.587.427.507.507.5010100.11
fazenda capoeirnha83.927.507.757.587.507.587.5010100.00
NA 82.587.427.427.507.337.257.6710100.00
various 83.087.837.337.507.677.587.6710100.10
NA 83.757.757.677.507.677.587.8310100.00
NA 82.177.587.677.427.177.177.7510100.08
NA 83.087.837.757.337.257.427.9210100.08
el papaturro 83.837.587.837.338.007.587.7510100.10
capoeirinha 83.927.587.677.677.757.587.8310100.11
In [44]:
dim( coffee_ratings_small )
  1. 636
  2. 11
In [45]:
coffee_ratings_small  %>% count(farm_name)
A tibble: 28 × 2
farm_namen
<fct><int>
agropecuaria quiagral 8
apollo estate 5
bethel 7
capoeirinha 10
cerro bueno 8
conquista / morito 11
doi tung development project 13
el morito 8
el papaturro 9
el sacramento 6
fazenda capoeirnha 13
finca medina 15
gamboa 7
kona pacific farmers cooperative 6
la esmeralda 7
la esperanza 7
la esperanza y anexos 5
la union monte verde 5
las cuchillas 6
las delicias 5
las merceditas 5
los hicaques 10
piamonte 5
rio verde 23
sethuraman estates 6
several 20
various 47
NA 359

Most of the farm names are missing.

Correlations:

In [46]:
options(repr.plot.width=15, repr.plot.height=15)
pairs( coffee_ratings_small  %>% select_if( is.numeric ) )
In [47]:
cor( coffee_ratings_small  %>% select_if( is.numeric ) )
A matrix: 10 × 10 of type dbl
total_cup_pointsaromaflavoraftertasteaciditybodybalanceuniformitysweetnessmoisture
total_cup_points 1.0000000 0.65681914 0.8015403 0.78574968 0.66177883 0.62615259 0.74319080.606135070.43420958-0.11188081
aroma 0.6568191 1.00000000 0.6810691 0.64683561 0.54205056 0.51869435 0.53700280.164897030.04948655-0.17164143
flavor 0.8015403 0.68106912 1.0000000 0.84272241 0.69323669 0.62554026 0.68230180.235414610.10927925-0.14365932
aftertaste 0.7857497 0.64683561 0.8427224 1.00000000 0.67870509 0.62304193 0.69491270.209659340.09521006-0.17606307
acidity 0.6617788 0.54205056 0.6932367 0.67870509 1.00000000 0.54324216 0.55206470.145101050.06173206-0.09009386
body 0.6261526 0.51869435 0.6255403 0.62304193 0.54324216 1.00000000 0.57754200.088953020.06447216-0.15556952
balance 0.7431908 0.53700281 0.6823018 0.69491270 0.55206465 0.57754195 1.00000000.255529790.10815129-0.19383086
uniformity 0.6061351 0.16489703 0.2354146 0.20965934 0.14510105 0.08895302 0.25552981.000000000.34773920 0.01767382
sweetness 0.4342096 0.04948655 0.1092792 0.09521006 0.06173206 0.06447216 0.10815130.347739201.00000000 0.09889003
moisture-0.1118808-0.17164143-0.1436593-0.17606307-0.09009386-0.15556952-0.19383090.017673820.09889003 1.00000000

Significant correlation among the predictors. Probably not all are necessary for a good model. Choose the following: uniformity, sweetness, moisture, flavor.

How to decide whether to use a model also with random slopes? Fit linear models for each subject!

In [48]:
options(repr.plot.width=15, repr.plot.height=3)
ggplot(coffee_ratings_small, aes(x=uniformity, y=total_cup_points, group=farm_name)) + 
    geom_smooth(method="lm", se=FALSE, linewidth=0.5)

`geom_smooth()` using formula = 'y ~ x'
In [49]:
options(repr.plot.width=15, repr.plot.height=3)
ggplot(coffee_ratings_small, aes(x=sweetness, y=total_cup_points, group=farm_name)) + 
    geom_smooth(method="lm", se=FALSE, linewidth=0.5)

`geom_smooth()` using formula = 'y ~ x'
In [50]:
options(repr.plot.width=15, repr.plot.height=3)
ggplot(coffee_ratings_small, aes(x=moisture, y=total_cup_points, group=farm_name)) + 
    geom_smooth(method="lm", se=FALSE, linewidth=0.5)

`geom_smooth()` using formula = 'y ~ x'
In [51]:
options(repr.plot.width=15, repr.plot.height=3)
ggplot(coffee_ratings_small, aes(x=flavor, y=total_cup_points, group=farm_name)) + 
    geom_smooth(method="lm", se=FALSE, linewidth=0.5)

`geom_smooth()` using formula = 'y ~ x'

Possibly, moisture and sweetness could be given individual slopes.

Fit and simulate two models¶

In [52]:
coffee_model_1 <- stan_glmer(
  total_cup_points ~ uniformity + sweetness + moisture + flavor + (1 | farm_name), 
  data = coffee_ratings_small, family = gaussian,
  prior_intercept = normal(250, 10, autoscale = TRUE),
  prior = normal(0, 1, autoscale = TRUE), 
  prior_aux = exponential(1, autoscale = TRUE),
  prior_covariance = decov(reg = 1, conc = 1, shape = 1, scale = 1),
  chains = 4, iter = 5000*2, seed = 84735
)
In [53]:
coffee_model_2 <- stan_glmer(
  total_cup_points ~ uniformity + sweetness + moisture + flavor + (moisture + sweetness | farm_name), 
  data = coffee_ratings_small, family = gaussian,
  prior_intercept = normal(250, 10, autoscale = TRUE),
  prior = normal(0, 1, autoscale = TRUE), 
  prior_aux = exponential(1, autoscale = TRUE),
  prior_covariance = decov(reg = 1, conc = 1, shape = 1, scale = 1),
  chains = 4, iter = 5000*2, seed = 84735
)

Warning message:
“There were 18 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.”
Warning message:
“Examine the pairs() plot to diagnose sampling problems
”

The second model takes longer on my laptop with four cores, if I wanted to include all predictors, this would take quite some time.. Regarding the warning: On https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup it says:

image.png

Let's check the diagnostics for both chains:

MCMC diagnostics¶

In [54]:
options(repr.plot.width=15, repr.plot.height=10)
mcmc_trace(coffee_model_1)
mcmc_dens_overlay(coffee_model_1)
options(repr.plot.width=25, repr.plot.height=20)
mcmc_acf(coffee_model_1)
neff_ratio(coffee_model_1)
rhat(coffee_model_1)
(Intercept)
0.96105
uniformity
0.9061
sweetness
0.81095
moisture
0.9627
flavor
0.6665
b[(Intercept) farm_name:agropecuaria_quiagral]
0.7088
b[(Intercept) farm_name:apollo_estate]
0.63165
b[(Intercept) farm_name:bethel]
0.90275
b[(Intercept) farm_name:capoeirinha]
1.14215
b[(Intercept) farm_name:cerro_bueno]
0.96395
b[(Intercept) farm_name:conquista_/_morito]
0.3612
b[(Intercept) farm_name:doi_tung_development_project]
0.87845
b[(Intercept) farm_name:el_morito]
1.15975
b[(Intercept) farm_name:el_papaturro]
1.11645
b[(Intercept) farm_name:el_sacramento]
0.9485
b[(Intercept) farm_name:fazenda_capoeirnha]
1.0215
b[(Intercept) farm_name:finca_medina]
0.86865
b[(Intercept) farm_name:gamboa]
0.5317
b[(Intercept) farm_name:kona_pacific_farmers_cooperative]
1.1473
b[(Intercept) farm_name:la_esmeralda]
1.2332
b[(Intercept) farm_name:la_esperanza]
1.2149
b[(Intercept) farm_name:la_esperanza_y_anexos]
1.24055
b[(Intercept) farm_name:la_union_monte_verde]
1.24415
b[(Intercept) farm_name:las_cuchillas]
0.81175
b[(Intercept) farm_name:las_delicias]
1.2122
b[(Intercept) farm_name:las_merceditas]
0.954
b[(Intercept) farm_name:los_hicaques]
1.0969
b[(Intercept) farm_name:piamonte]
1.2971
b[(Intercept) farm_name:rio_verde]
0.53975
b[(Intercept) farm_name:sethuraman_estates]
0.62595
b[(Intercept) farm_name:several]
0.9423
b[(Intercept) farm_name:various]
0.7126
sigma
0.7494
Sigma[farm_name:(Intercept),(Intercept)]
0.3188
(Intercept)
1.0001447719945
uniformity
1.00036696897546
sweetness
1.00070041410262
moisture
0.999998074523331
flavor
1.00004094789613
b[(Intercept) farm_name:agropecuaria_quiagral]
1.00000353590096
b[(Intercept) farm_name:apollo_estate]
0.999957929246099
b[(Intercept) farm_name:bethel]
0.999977285065357
b[(Intercept) farm_name:capoeirinha]
0.999944994791927
b[(Intercept) farm_name:cerro_bueno]
1.00009273455109
b[(Intercept) farm_name:conquista_/_morito]
1.00006720788921
b[(Intercept) farm_name:doi_tung_development_project]
0.999903778115077
b[(Intercept) farm_name:el_morito]
0.999944900907315
b[(Intercept) farm_name:el_papaturro]
1.00006005649433
b[(Intercept) farm_name:el_sacramento]
1.00000911660307
b[(Intercept) farm_name:fazenda_capoeirnha]
0.999969380990973
b[(Intercept) farm_name:finca_medina]
1.00024127601169
b[(Intercept) farm_name:gamboa]
0.99987430745209
b[(Intercept) farm_name:kona_pacific_farmers_cooperative]
0.999897750250558
b[(Intercept) farm_name:la_esmeralda]
0.999925996940665
b[(Intercept) farm_name:la_esperanza]
0.999932256746219
b[(Intercept) farm_name:la_esperanza_y_anexos]
0.99991843494923
b[(Intercept) farm_name:la_union_monte_verde]
0.999956898202949
b[(Intercept) farm_name:las_cuchillas]
1.00001035254489
b[(Intercept) farm_name:las_delicias]
0.999857790684174
b[(Intercept) farm_name:las_merceditas]
0.999902807615587
b[(Intercept) farm_name:los_hicaques]
0.999916502974856
b[(Intercept) farm_name:piamonte]
0.999956408185194
b[(Intercept) farm_name:rio_verde]
1.00007254677106
b[(Intercept) farm_name:sethuraman_estates]
1.00064542378429
b[(Intercept) farm_name:several]
1.00010536912856
b[(Intercept) farm_name:various]
1.00002000340576
sigma
0.99993283478065
Sigma[farm_name:(Intercept),(Intercept)]
1.00030618584635

Diagnostics look good.

In [55]:
options(repr.plot.width=15, repr.plot.height=25)
mcmc_trace(coffee_model_2)
mcmc_dens_overlay(coffee_model_2)
mcmc_acf(coffee_model_2)
neff_ratio(coffee_model_2)
rhat(coffee_model_2)
(Intercept)
1.08505
uniformity
1.18815
sweetness
1.05285
moisture
1.2604
flavor
0.93525
b[(Intercept) farm_name:agropecuaria_quiagral]
0.5266
b[moisture farm_name:agropecuaria_quiagral]
0.53495
b[sweetness farm_name:agropecuaria_quiagral]
0.7959
b[(Intercept) farm_name:apollo_estate]
0.44415
b[moisture farm_name:apollo_estate]
0.41435
b[sweetness farm_name:apollo_estate]
0.7474
b[(Intercept) farm_name:bethel]
0.6638
b[moisture farm_name:bethel]
0.60895
b[sweetness farm_name:bethel]
1.02055
b[(Intercept) farm_name:capoeirinha]
0.93845
b[moisture farm_name:capoeirinha]
0.8253
b[sweetness farm_name:capoeirinha]
1.32635
b[(Intercept) farm_name:cerro_bueno]
0.71295
b[moisture farm_name:cerro_bueno]
0.68925
b[sweetness farm_name:cerro_bueno]
1.1131
b[(Intercept) farm_name:conquista_/_morito]
0.29565
b[moisture farm_name:conquista_/_morito]
0.30045
b[sweetness farm_name:conquista_/_morito]
0.4942
b[(Intercept) farm_name:doi_tung_development_project]
0.7176
b[moisture farm_name:doi_tung_development_project]
0.61345
b[sweetness farm_name:doi_tung_development_project]
1.0349
b[(Intercept) farm_name:el_morito]
0.90925
b[moisture farm_name:el_morito]
0.80495
b[sweetness farm_name:el_morito]
1.17875
b[(Intercept) farm_name:el_papaturro]
1.05835
b[moisture farm_name:el_papaturro]
0.8508
b[sweetness farm_name:el_papaturro]
1.2121
b[(Intercept) farm_name:el_sacramento]
0.79275
b[moisture farm_name:el_sacramento]
0.6977
b[sweetness farm_name:el_sacramento]
1.09965
b[(Intercept) farm_name:fazenda_capoeirnha]
0.85495
b[moisture farm_name:fazenda_capoeirnha]
0.8731
b[sweetness farm_name:fazenda_capoeirnha]
1.17715
b[(Intercept) farm_name:finca_medina]
0.6888
b[moisture farm_name:finca_medina]
0.7233
b[sweetness farm_name:finca_medina]
0.962
b[(Intercept) farm_name:gamboa]
0.39135
b[moisture farm_name:gamboa]
0.36645
b[sweetness farm_name:gamboa]
0.57275
b[(Intercept) farm_name:kona_pacific_farmers_cooperative]
0.9374
b[moisture farm_name:kona_pacific_farmers_cooperative]
0.8832
b[sweetness farm_name:kona_pacific_farmers_cooperative]
1.2761
b[(Intercept) farm_name:la_esmeralda]
0.97815
b[moisture farm_name:la_esmeralda]
0.8361
b[sweetness farm_name:la_esmeralda]
1.2681
b[(Intercept) farm_name:la_esperanza]
0.7806
b[moisture farm_name:la_esperanza]
0.8406
b[sweetness farm_name:la_esperanza]
1.24595
b[(Intercept) farm_name:la_esperanza_y_anexos]
0.8026
b[moisture farm_name:la_esperanza_y_anexos]
0.802
b[sweetness farm_name:la_esperanza_y_anexos]
1.11335
b[(Intercept) farm_name:la_union_monte_verde]
0.88595
b[moisture farm_name:la_union_monte_verde]
0.8622
b[sweetness farm_name:la_union_monte_verde]
1.2965
b[(Intercept) farm_name:las_cuchillas]
0.5789
b[moisture farm_name:las_cuchillas]
0.54485
b[sweetness farm_name:las_cuchillas]
0.92865
b[(Intercept) farm_name:las_delicias]
1.0194
b[moisture farm_name:las_delicias]
0.8341
b[sweetness farm_name:las_delicias]
1.21465
b[(Intercept) farm_name:las_merceditas]
0.7979
b[moisture farm_name:las_merceditas]
0.62775
b[sweetness farm_name:las_merceditas]
1.08125
b[(Intercept) farm_name:los_hicaques]
0.8622
b[moisture farm_name:los_hicaques]
0.86195
b[sweetness farm_name:los_hicaques]
1.206
b[(Intercept) farm_name:piamonte]
0.8999
b[moisture farm_name:piamonte]
0.85695
b[sweetness farm_name:piamonte]
1.25015
b[(Intercept) farm_name:rio_verde]
0.3926
b[moisture farm_name:rio_verde]
0.35455
b[sweetness farm_name:rio_verde]
0.60675
b[(Intercept) farm_name:sethuraman_estates]
0.56345
b[moisture farm_name:sethuraman_estates]
0.5369
b[sweetness farm_name:sethuraman_estates]
0.85025
b[(Intercept) farm_name:several]
0.9247
b[moisture farm_name:several]
0.8257
b[sweetness farm_name:several]
1.0299
b[(Intercept) farm_name:various]
0.67325
b[moisture farm_name:various]
0.7865
b[sweetness farm_name:various]
0.8879
sigma
0.88665
Sigma[farm_name:(Intercept),(Intercept)]
0.6373
Sigma[farm_name:moisture,(Intercept)]
1.07195
Sigma[farm_name:sweetness,(Intercept)]
0.5102
Sigma[farm_name:moisture,moisture]
0.48565
Sigma[farm_name:sweetness,moisture]
0.23595
Sigma[farm_name:sweetness,sweetness]
0.3677
(Intercept)
0.99988293226406
uniformity
0.999994399715048
sweetness
1.00006165018381
moisture
0.999894559777911
flavor
1.00001076158283
b[(Intercept) farm_name:agropecuaria_quiagral]
1.00029957325701
b[moisture farm_name:agropecuaria_quiagral]
1.00010798333829
b[sweetness farm_name:agropecuaria_quiagral]
0.999966902183615
b[(Intercept) farm_name:apollo_estate]
0.999884295535258
b[moisture farm_name:apollo_estate]
0.99994364253767
b[sweetness farm_name:apollo_estate]
0.999991226126833
b[(Intercept) farm_name:bethel]
0.999954327880675
b[moisture farm_name:bethel]
0.999960711430095
b[sweetness farm_name:bethel]
1.00005040949724
b[(Intercept) farm_name:capoeirinha]
0.999898197493059
b[moisture farm_name:capoeirinha]
0.999936064313419
b[sweetness farm_name:capoeirinha]
0.99986517056971
b[(Intercept) farm_name:cerro_bueno]
1.0000966633813
b[moisture farm_name:cerro_bueno]
1.0001957892794
b[sweetness farm_name:cerro_bueno]
1.00000752135181
b[(Intercept) farm_name:conquista_/_morito]
1.00017754665898
b[moisture farm_name:conquista_/_morito]
1.00014616194474
b[sweetness farm_name:conquista_/_morito]
1.00009310407862
b[(Intercept) farm_name:doi_tung_development_project]
1.00057147020839
b[moisture farm_name:doi_tung_development_project]
1.00001472961773
b[sweetness farm_name:doi_tung_development_project]
1.00012562836054
b[(Intercept) farm_name:el_morito]
1.00033105067131
b[moisture farm_name:el_morito]
1.00012588627906
b[sweetness farm_name:el_morito]
0.999998496561761
b[(Intercept) farm_name:el_papaturro]
0.999977386163742
b[moisture farm_name:el_papaturro]
0.999915236267374
b[sweetness farm_name:el_papaturro]
1.00008497205081
b[(Intercept) farm_name:el_sacramento]
0.999884330111701
b[moisture farm_name:el_sacramento]
1.00038281750128
b[sweetness farm_name:el_sacramento]
0.999979866391964
b[(Intercept) farm_name:fazenda_capoeirnha]
1.00016913416359
b[moisture farm_name:fazenda_capoeirnha]
0.999936892920424
b[sweetness farm_name:fazenda_capoeirnha]
0.999903749363805
b[(Intercept) farm_name:finca_medina]
0.999896233405066
b[moisture farm_name:finca_medina]
1.00006059361785
b[sweetness farm_name:finca_medina]
0.999996621553092
b[(Intercept) farm_name:gamboa]
0.999960321351392
b[moisture farm_name:gamboa]
1.00009950515386
b[sweetness farm_name:gamboa]
0.99998712462226
b[(Intercept) farm_name:kona_pacific_farmers_cooperative]
1.00029809385987
b[moisture farm_name:kona_pacific_farmers_cooperative]
0.999973304936992
b[sweetness farm_name:kona_pacific_farmers_cooperative]
0.999914115809875
b[(Intercept) farm_name:la_esmeralda]
1.00010847003517
b[moisture farm_name:la_esmeralda]
1.00035073078166
b[sweetness farm_name:la_esmeralda]
0.999988086662427
b[(Intercept) farm_name:la_esperanza]
0.999899544076697
b[moisture farm_name:la_esperanza]
0.999889325741995
b[sweetness farm_name:la_esperanza]
0.999838284560986
b[(Intercept) farm_name:la_esperanza_y_anexos]
0.999914388647449
b[moisture farm_name:la_esperanza_y_anexos]
0.999921669986585
b[sweetness farm_name:la_esperanza_y_anexos]
0.999886670167431
b[(Intercept) farm_name:la_union_monte_verde]
1.00006129679379
b[moisture farm_name:la_union_monte_verde]
0.999932199005353
b[sweetness farm_name:la_union_monte_verde]
0.999893138470969
b[(Intercept) farm_name:las_cuchillas]
1.00041317585073
b[moisture farm_name:las_cuchillas]
1.00035781871895
b[sweetness farm_name:las_cuchillas]
1.00005507893655
b[(Intercept) farm_name:las_delicias]
1.00010290752743
b[moisture farm_name:las_delicias]
1.0001157694926
b[sweetness farm_name:las_delicias]
0.999976634924647
b[(Intercept) farm_name:las_merceditas]
0.999941386271519
b[moisture farm_name:las_merceditas]
1.0000707045289
b[sweetness farm_name:las_merceditas]
0.999985328329294
b[(Intercept) farm_name:los_hicaques]
1.00001928948171
b[moisture farm_name:los_hicaques]
1.00026721142856
b[sweetness farm_name:los_hicaques]
0.999886564060915
b[(Intercept) farm_name:piamonte]
0.999995138241298
b[moisture farm_name:piamonte]
1.00000987659665
b[sweetness farm_name:piamonte]
1.00005393290013
b[(Intercept) farm_name:rio_verde]
0.999925530197326
b[moisture farm_name:rio_verde]
1.00021132292098
b[sweetness farm_name:rio_verde]
0.999982542084742
b[(Intercept) farm_name:sethuraman_estates]
1.00065244782644
b[moisture farm_name:sethuraman_estates]
1.00011918623267
b[sweetness farm_name:sethuraman_estates]
1.00003482901908
b[(Intercept) farm_name:several]
1.00014441887811
b[moisture farm_name:several]
1.00033024748864
b[sweetness farm_name:several]
1.00003835843767
b[(Intercept) farm_name:various]
1.00016601748738
b[moisture farm_name:various]
1.00003591567164
b[sweetness farm_name:various]
1.00010671631935
sigma
0.99994294408406
Sigma[farm_name:(Intercept),(Intercept)]
1.00011102023504
Sigma[farm_name:moisture,(Intercept)]
0.999932742974669
Sigma[farm_name:sweetness,(Intercept)]
1.00008435340285
Sigma[farm_name:moisture,moisture]
1.00032999292601
Sigma[farm_name:sweetness,moisture]
1.00047839807522
Sigma[farm_name:sweetness,sweetness]
0.999894630140094

There are some spikes (probably the divergences mentioned in the warning) in the trace plots. However the density plots and the rhat values seem to look fine. To do this properly, one would need to tweak the MCMC here, however I ignore this and continue.

Posterior predictive check¶

In [56]:
options(repr.plot.width=15, repr.plot.height=5)
pp_check(coffee_model_1)
In [57]:
options(repr.plot.width=15, repr.plot.height=5)
pp_check(coffee_model_2)

Both models look more or less ok, the actual data distribution is a bit more narrow than the normal distribution.

Compare model performances and select best model¶

In [58]:
prediction_summary(model = coffee_model_1, data = coffee_ratings_small %>% na.omit())
A data.frame: 1 × 4
maemae_scaledwithin_50within_95
<dbl><dbl><dbl><dbl>
0.4928290.49018730.60288810.967509
In [59]:
prediction_summary(model = coffee_model_2, data = coffee_ratings_small %>% na.omit())
A data.frame: 1 × 4
maemae_scaledwithin_50within_95
<dbl><dbl><dbl><dbl>
0.48720590.48390240.60649820.967509

ELPD:

In [60]:
elpd1 <- loo(coffee_model_1, k_threshold=0.7)
elpd2 <- loo(coffee_model_2, k_threshold=0.7)
loo_compare(elpd1, elpd2)

1 problematic observation(s) found.
Model will be refit 1 times.


Fitting model 1 out of 1 (leaving out observation 213)

1 problematic observation(s) found.
Model will be refit 1 times.


Fitting model 1 out of 1 (leaving out observation 213)
A compare.loo: 2 × 8 of type dbl
elpd_diffse_diffelpd_loose_elpd_loop_loose_p_loolooicse_looic
coffee_model_2 0.0000000000.0000000-399.164822.9940024.198246.241245798.329545.98801
coffee_model_1-0.0092634330.3397638-399.174022.9229623.898636.114573798.348145.84593

No significant difference between the models, only use the simpler model 1 from now on.

Global model parameters¶

In [61]:
tidy(coffee_model_1, effects = "fixed", conf.int = TRUE, conf.level = 0.80)
A tibble: 5 × 5
termestimatestd.errorconf.lowconf.high
<chr><dbl><dbl><dbl><dbl>
(Intercept)13.8236471.859174411.384129916.17473545
uniformity 1.8725900.1208561 1.7170889 2.02891949
sweetness 1.0913690.1381829 0.9116501 1.26857769
moisture -1.7323041.4074539-3.5386459 0.08841798
flavor 5.2188790.1932896 4.9743177 5.47074034

Global average model:

Total cup points = 13.8 + uniformity * 1.9 + sweetness * 1.09 - 1.73 moisture + 5.2 * flavor

However the negative relationship with moisture seems not to be significant at the 80% level. The strongest contribution comes from flavor: Every point in flavor gives 5 points in total cup points.

Variability¶

Which variability is larger, between-group (between all the farms) or within-group? (among one farm) ?

In [62]:
tidy(coffee_model_1, effects = "ran_pars")
A tibble: 2 × 3
termgroupestimate
<chr><chr><dbl>
sd_(Intercept).farm_namefarm_name0.3490031
sd_Observation.Residual Residual 0.9753836
In [63]:
sigma_y <- 0.97
sigma_0 <- 0.35

Proportion of variability explained by within-group variability:

In [64]:
sigma_y^2/(sigma_y^2 + sigma_0^2)
0.884803460598082

Proportion of variability explained by between-group variability:

In [65]:
sigma_0^2/(sigma_y^2 + sigma_0^2)
0.115196539401918

It appears that the coffee farm alone is not really an indicator for good coffee.

In [66]:
ggplot( coffee_ratings_small, aes(x=flavor, y=total_cup_points, fill=farm_name) ) + geom_boxplot()

Which is the best farm, which the worst?¶

As explained above, the within-group variability dominates over between-group variability, so there are only small differences. Apparently stan_glmer() knows how to deal with missing data (most of the farm names are missing), I ran both simulations with and without na.omit() and got the same results.

In [67]:
coffee_summary <- coffee_model_1 %>%
  spread_draws(`(Intercept)`, b[,farm_name]) %>% # extract dataframe
  mutate(farm_name_intercept = `(Intercept)` + b) %>% # compute individual runner intercepts
  select(-`(Intercept)`, -b) %>% # drop global intercept and b
  median_qi(.width = 0.80) %>%  # compute credible intervals
  select(farm_name, farm_name_intercept, .lower, .upper)

coffee_summary %>% arrange( desc(farm_name_intercept) )
A tibble: 27 × 4
farm_namefarm_name_intercept.lower.upper
<chr><dbl><dbl><dbl>
farm_name:conquista_/_morito 14.4220512.0925516.67708
farm_name:gamboa 14.2674511.8637216.60838
farm_name:apollo_estate 14.1577811.7493616.51272
farm_name:sethuraman_estates 14.0943011.7513316.39570
farm_name:doi_tung_development_project 14.0249811.5766916.36636
farm_name:los_hicaques 13.8926911.4254016.26750
farm_name:capoeirinha 13.8647311.3926816.25867
farm_name:el_morito 13.8621711.4009116.23128
farm_name:la_union_monte_verde 13.8548911.4012916.24049
farm_name:fazenda_capoeirnha 13.8424411.3591016.20933
farm_name:la_esmeralda 13.8112711.3452916.20064
farm_name:las_delicias 13.8028711.3328616.18417
farm_name:piamonte 13.8017211.3372516.20597
farm_name:kona_pacific_farmers_cooperative13.7948911.2967416.19857
farm_name:various 13.7836411.3768716.13317
farm_name:el_papaturro 13.7836111.2937816.19005
farm_name:la_esperanza 13.7799011.3181916.17485
farm_name:several 13.7637511.2921816.15315
farm_name:la_esperanza_y_anexos 13.7343611.2303016.15126
farm_name:finca_medina 13.6869211.1911516.10384
farm_name:el_sacramento 13.6811511.1891016.10473
farm_name:las_merceditas 13.6680711.1843416.07224
farm_name:cerro_bueno 13.6199111.2090915.98864
farm_name:bethel 13.6100311.1160216.02990
farm_name:las_cuchillas 13.5998011.0860016.02897
farm_name:agropecuaria_quiagral 13.5472911.0521715.96046
farm_name:rio_verde 13.4411310.9519015.85587

Exercise 17.15¶

Prior elicitation¶

For the chosen hierarchical random intercept model:

\begin{eqnarray*} Y_{ij}|\beta_{0j},\beta_1,\sigma_y &\sim& N(\mu_{ij}, \sigma_y^2) \quad \text{with} \quad \mu_{ij} = \beta_{0j} + \beta_{1} X _{ij}\\ \beta_{0j} |\beta_0, \sigma_0 &\sim & N\left( \beta_0, \sigma_0^2 \right)\\ \beta_{0c} &=& N(m_0, s_0^2) \\ \beta_1 &=& N(0, s_1^2)\\ \sigma_y &=& \text{Exp}(l_y) \\ \sigma_0 &=& \text{Exp}(l_0) \end{eqnarray*}

Average reaction time is probably less than half a second for all the participants, however also larger than 100 ms:

In [68]:
plot_normal( mean=250, sd=100 )

With my very weak knowledge: The increase of reaction time per day is certainly only a part of the full baseline reaction time and probably not more than 50$\mu$s per day on average. A decrease in reaction time should be an expection. This is difficult to model, as choosing too large slopes here might lead to many prior models with a negative intercept. Do do this analysis in a theoretically more sound way, one should use a prior distribution that makes sure that intercepts are only positive, however this normal model appears to be useful nevertheless.

In [69]:
plot_normal( mean=25, sd=10 )

As mean within-group variability, I expect around $\sigma_y \sim 100\mu$s: $l_y = 1/100$. As mean between-group variability I expect around $\sigma_0 \sim 200\mu$s: $l_0 = 1/200$. Now things get a bit fuzzier and I am not perfectly sure if this is the right solution. If I understood things correctly, the reg and conc parameters of decov are not important here, since we do not have to model $\sigma_1$ and hence not a multivariate normal distribution, but only the boundary case of an univariate one. The shape and scale parameters are related to the prior for $\sigma_0$ (that should be equal to $\tau$): $\text{Gamma}(1,\lambda) = \text{Exp}(\lambda)$. Thus I use decov(reg=1, conc=1, shape=1, scale=1/200).

Prior simulation¶

In [70]:
sleep_model_priors <- stan_glmer(
  Reaction ~ Days + (1 | Subject), 
  data = sleepstudy, family = gaussian,
  prior_intercept = normal(250, 100),
  prior = normal(25, 10), 
  prior_aux = exponential(1/100),
  prior_covariance = decov(reg = 1, conc = 1, shape = 1, scale = 1/200),
  chains = 4, iter = 5000*2, seed = 84735,
  prior_PD = TRUE
)

Checkout prior models¶

In [71]:
sleepstudy %>%
  add_fitted_draws(sleep_model_priors, n = 200) %>%
  ggplot(aes(x = Days, y = Reaction)) +
    geom_line(aes(y = .value, group = .draw), alpha = 0.15)

Warning message:
“`fitted_draws` and `add_fitted_draws` are deprecated as their names were confusing.
- Use [add_]epred_draws() to get the expectation of the posterior predictive.
- Use [add_]linpred_draws() to get the distribution of the linear predictor.
- For example, you used [add_]fitted_draws(..., scale = "response"), which
  means you most likely want [add_]epred_draws(...).
NOTE: When updating to the new functions, note that the `model` parameter is now
  named `object` and the `n` parameter is now named `ndraws`.”

Simulate a few datasets:

In [72]:
sleepstudy %>%
  add_predicted_draws(sleep_model_priors, n = 12) %>%
  ggplot(aes(x = Days, y = Reaction)) +
    geom_point(aes(y = .prediction, group = .draw)) + 
    facet_wrap(~ .draw)

Warning message:
“
In add_predicted_draws(): The `n` argument is a deprecated alias for `ndraws`.
Use the `ndraws` argument instead.
See help("tidybayes-deprecated").
”

Looks more or less reasonable.